###############################################
# This code produces 
# 1) Correlation tests 
# 2) Figures 1 ~ 7
# 3) Tables 1 ~ 5
# 4) Effect sizes 
###############################################

library(tidyverse)
library(bucky)
library(ggeffects)
library(ggpubr)
library(broom)
library(stargazer)

setwd() # set working directory to the folder with the replication materials

# Create list of countries with radical right and radical left parties
partylist <- read.csv("party_list.csv")
c_rr <- partylist %>%
  filter(has_seats==1, populist == 1, far_right == 1) %>%
  dplyr::select(cntry) %>% distinct() %>% pull() %>%
  c(.,"United Kingdom") # Add UK (UKIP)

c_rl <- partylist %>%
  filter(has_seats==1, populist == 1, far_left == 1) %>%
  filter(cntry != "United Kingdom") %>% # it's incorrectly included because of Sinn Féin
  dplyr::select(cntry) %>% distinct() %>% pull()

# Load cleaned ESS data
dat <- read.csv("ESS_cleaned.csv") 

# Set theme for figures
theme_set(theme_bw(base_size = 16) +
            theme(axis.line = element_line(color='black'),
                  plot.background = element_blank(),
                  panel.grid.minor = element_blank(),
                  panel.grid.major = element_blank()))

# Correlation ---------------------------------------------------------------
# income decile and unfairness of income of their own income
cor.test(dat$income, dat$inc_unfair,
         method = "pearson", use = "complete.obs")

# income decile and unfairness of income of the top 10%
cor.test(dat$income, dat$top_inc_unfair,
         method = "pearson", use = "complete.obs")

# income decile and unfairness of income of the bottom 10%
cor.test(dat$income, dat$bttm_inc_unfair,
         method = "pearson", use = "complete.obs")

# Tables ---------------------------------------------------------------

##########################
# Table 1
##########################

# Column 1
t1_1 <- lm( inc_unfair ~ RTI + netjobreduction_manu_scaled + netjobreduction_non_manu_scaled + 
              college + female + age + income + union + unemployed + 
              citizen + religion + urban + 
              college_share + male_share + unemp_rate + imm_share + cntry,
            data = dat, weights = dweight) 
t1_1 <- robustify(t1_1 , cluster = code_2003)

# Column 2
t1_2 <- lm( top_inc_unfair ~ RTI + netjobreduction_manu_scaled + netjobreduction_non_manu_scaled + 
              college + female + age + income + union + unemployed + 
              citizen + religion + urban + 
              college_share + male_share + unemp_rate + imm_share + cntry,
            data = dat, weights = dweight) 
t1_2 <- robustify(t1_2, cluster = code_2003)

# Column 3
t1_3 <- lm( bttm_inc_unfair ~ RTI + netjobreduction_manu_scaled + netjobreduction_non_manu_scaled + 
              college + female + age + income + union + unemployed + 
              citizen + religion + urban + 
              college_share + male_share + unemp_rate + imm_share + cntry,
            data = dat, weights = dweight) 
t1_3 <- robustify(t1_3, cluster = code_2003)

# Column 4
t1_4 <- lm( job_self_unfair ~RTI + netjobreduction_manu_scaled + netjobreduction_non_manu_scaled + 
              college + female + age + income + union + unemployed + 
              citizen + religion + urban + 
              college_share + male_share + unemp_rate + imm_share + cntry,
            data = dat, weights = dweight) 
t1_4 <- robustify(t1_4, cluster = code_2003)

# Column 5
t1_5 <- lm( job_other_unfair ~ RTI + netjobreduction_manu_scaled + netjobreduction_non_manu_scaled + 
              college + female + age + income + union + unemployed + 
              citizen + religion + urban + 
              college_share + male_share + unemp_rate + imm_share + cntry,
            data = dat, weights = dweight) 
t1_5 <- robustify(t1_5, cluster = code_2003)

stargazer(t1_1,t1_2,t1_3,t1_4,t1_5, type = "text")

##########################
# Table 2
##########################

# Column 1
t2_1 <- lm( redist ~ inc_unfair + top_inc_unfair + bttm_inc_unfair + 
              college + female + age + income + union + unemployed + 
              citizen + religion + urban + cntry,
            data = dat, weights = dweight) 
t2_1 <- robustify(t2_1 , cluster = cntry)

# Column 2 
t2_2 <- lm( redist ~ job_self_unfair + job_other_unfair +
              college + female + age + income + union + unemployed + 
              citizen + religion + urban + cntry,
            data = dat, weights = dweight) 
t2_2 <- robustify(t2_2 , cluster = cntry)

# Column 3
t2_3 <- lm( imm_culture ~ inc_unfair + top_inc_unfair + bttm_inc_unfair + 
              college + female + age + income + union + unemployed + 
              citizen + religion + urban + cntry,
            data = dat, weights = dweight) 
t2_3 <- robustify(t2_3 , cluster = cntry)

# Column 4
t2_4 <- lm( imm_culture ~ job_self_unfair + job_other_unfair +
              college + female + age + income + union + unemployed + 
              citizen + religion + urban + cntry,
            data = dat, weights = dweight) 
t2_4 <- robustify(t2_4 , cluster = cntry)

# Combined
stargazer(t2_1,t2_2,t2_3,t2_4, type = "text")

##########################
# Table 3
##########################

# Column 1
t3_1 <- glm( far_left ~ inc_unfair + top_inc_unfair + bttm_inc_unfair + 
               college + female + age + income + union + unemployed + 
               citizen + religion + urban + cntry,
             data = dat  %>% filter(cntry %in% c_rl),
             weights = dweight, family=binomial(link=logit)) 
t3_1 <- robustify(t3_1 , cluster = cntry)

# Column 2
t3_2 <- glm( far_left ~ job_self_unfair + job_other_unfair +
               college + female + age + income + union + unemployed + 
               citizen + religion + urban + cntry,
             data = dat %>% filter(cntry %in% c_rl),
             weights = dweight, family=binomial(link=logit)) 
t3_2 <- robustify(t3_2 , cluster = cntry)

# Column 3
t3_3 <- glm( far_right ~ inc_unfair + top_inc_unfair + bttm_inc_unfair + 
               college + female + age + income + union + unemployed + 
               citizen + religion + urban + cntry,
             data = dat %>% filter(cntry %in% c(c_rr, "United Kingdom")),
             weights = dweight, family=binomial(link=logit)) 
t3_3 <- robustify(t3_3 , cluster = cntry)

# Column 4 
t3_4 <- glm( far_right ~ job_self_unfair + job_other_unfair +
               college + female + age + income + union + unemployed + 
               citizen + religion + urban + cntry,
             data = dat %>% filter(cntry %in% c(c_rr)),
             weights = dweight, family=binomial(link=logit)) 
t3_4 <- robustify(t3_4 , cluster = cntry)

# Combined
stargazer(t3_1,t3_2,t3_3,t3_4, type = "text")

##########################
# Table 4
##########################

# Column 1
t4_1 <- glm( center_left ~ inc_unfair + top_inc_unfair + bttm_inc_unfair + 
               college + female + age + income + union + unemployed + 
               citizen + religion + urban + cntry,
             data = dat,
             weights = dweight, family=binomial(link=logit)) 
t4_1 <- robustify(t4_1 , cluster = cntry)

# Column 2
t4_2 <- glm( center_left ~ job_self_unfair + job_other_unfair +
               college + female + age + income + union + unemployed + 
               citizen + religion + urban + cntry,
             data = dat,
             weights = dweight, family=binomial(link=logit)) 
t4_2 <- robustify(t4_2 , cluster = cntry)

# Column 3
t4_3 <- glm( center_right ~ inc_unfair + top_inc_unfair + bttm_inc_unfair + 
               college + female + age + income + union + unemployed + 
               citizen + religion + urban + cntry,
             data = dat ,
             weights = dweight, family=binomial(link=logit)) 
t4_3 <- robustify(t4_3 , cluster = cntry)

# Column 4 
t4_4 <- glm( center_right ~ job_self_unfair + job_other_unfair +
               college + female + age + income + union + unemployed + 
               citizen + religion + urban + cntry,
             data = dat ,
             weights = dweight, family=binomial(link=logit)) 
t4_4 <- robustify(t4_4 , cluster = cntry)

# Column 5
t4_5 <- glm( voted ~ inc_unfair + top_inc_unfair + bttm_inc_unfair + 
               college + female + age + income + union + unemployed + 
               citizen + religion + urban + cntry,
             data = dat,
             weights = dweight, family=binomial(link=logit)) 
t4_5 <- robustify(t4_5, cluster = cntry)

# Column 6
t4_6 <- glm( voted ~ job_self_unfair + job_other_unfair +
               college + female + age + income + union + unemployed + 
               citizen + religion + urban + cntry,
             data = dat,
             weights = dweight, family=binomial(link=logit)) 
t4_6 <- robustify(t4_6, cluster = cntry)

# Combined
stargazer(t4_1,t4_2,t4_3,t4_4,t4_5,t4_6, type = "text")

##########################
# Table 5
##########################

# Standardize variables for comparison
dat <- dat %>%
  mutate(abs_deprivation_std = (abs_deprivation - weighted.mean(abs_deprivation, w=dweight,na.rm=T))/sqrt(wtd.var(abs_deprivation, w = dweight, na.rm=T)),
         authoritarian_std = (authoritarian - weighted.mean(authoritarian, w=dweight,na.rm=T))/sqrt(wtd.var(authoritarian, w = dweight, na.rm=T)),
         inc_unfair_std = (inc_unfair - weighted.mean(inc_unfair, w=dweight,na.rm=T))/sqrt(wtd.var(inc_unfair, w = dweight, na.rm=T)),
         top_inc_unfair_std = (top_inc_unfair - weighted.mean(top_inc_unfair, w=dweight,na.rm=T))/sqrt(wtd.var(top_inc_unfair, w = dweight, na.rm=T)),
         bttm_inc_unfair_std = (bttm_inc_unfair - weighted.mean(bttm_inc_unfair, w=dweight,na.rm=T))/sqrt(wtd.var(bttm_inc_unfair, w = dweight, na.rm=T)),
         job_self_unfair_std = (job_self_unfair - weighted.mean(job_self_unfair, w=dweight,na.rm=T))/sqrt(wtd.var(job_self_unfair, w = dweight, na.rm=T)),
         job_other_unfair_std = (job_other_unfair - weighted.mean(job_other_unfair, w=dweight,na.rm=T))/sqrt(wtd.var(job_other_unfair, w = dweight, na.rm=T))
  )

# Column 1
t5_1 <- glm( far_left ~ inc_unfair_std + top_inc_unfair_std + bttm_inc_unfair_std + 
               abs_deprivation_std + authoritarian_std +
               college + female + age + income + union + unemployed + 
               citizen + religion + urban + cntry,
             data = dat  %>% filter(cntry %in% c_rl),
             weights = dweight, family=binomial(link=logit)) 
t5_1 <- robustify(t5_1, cluster = cntry)

# Column 2
t5_2 <- glm( far_left ~ job_self_unfair_std + job_other_unfair_std + 
               abs_deprivation_std +  authoritarian_std  +
               college + female + age + income + union + unemployed + 
               citizen + religion + urban + cntry,
             data = dat %>% filter(cntry %in% c_rl),
             weights = dweight, family=binomial(link=logit)) 
t5_2 <- robustify(t5_2 , cluster = cntry)

# Column 3
t5_3 <- glm( far_right ~ inc_unfair_std + top_inc_unfair_std + bttm_inc_unfair_std + 
               abs_deprivation_std +  authoritarian_std  +
               college + female + age + income + union + unemployed + 
               citizen + religion + urban + cntry,
             data = dat %>% filter(cntry %in% c(c_rr)),
             weights = dweight, family=binomial(link=logit)) 
t5_3 <- robustify(t5_3, cluster = cntry)

# Column 4 
t5_4 <- glm( far_right ~ job_self_unfair_std + job_other_unfair_std +
               abs_deprivation_std +  authoritarian_std  +
               college + female + age + income + union + unemployed + 
               citizen + religion + urban + cntry,
             data = dat %>% filter(cntry %in% c(c_rr)),
             weights = dweight, family=binomial(link=logit)) 
t5_4 <- robustify(t5_4, cluster = cntry)

stargazer(t5_1,t5_2,t5_3,t5_4, type = "text")

# Figures ---------------------------------------------------------------

## Create 'figure' folder to save plots

##########################
# Figure 1
##########################
(f1 <- dat %>%
   filter(is.na(income) == 0, is.na(inc_unfair) == 0) %>%
   ggplot(aes(x=inc_unfair)) +
   geom_bar(aes(stat="identity")) +
   facet_wrap(~income, ncol = 5) + 
   scale_x_continuous(breaks=-4:4) +
   xlab("Perceived Unfairness of Income for Self") + ylab("Count"))

ggsave("figure/Figure1.jpeg", f1, width = 10, height = 5)


# Appendix Figures B1-B4
(b1 <- dat %>%
   filter(is.na(income) == 0, is.na(bttm_inc_unfair) == 0) %>%
   ggplot(aes(x=bttm_inc_unfair)) +
   geom_bar(aes(stat="identity")) +
   facet_wrap(~income, ncol = 5) + 
   scale_x_continuous(breaks=-4:4) +
   xlab("Perceived Unfairness of Bottom 10% Income") + ylab("Count"))

ggsave("figure/FigureB1_Appendix.jpeg", b1, width = 10, height = 5)

# Figure B2
(b2 <- dat %>%
    filter(is.na(income) == 0, is.na(top_inc_unfair) == 0) %>%
    ggplot(aes(x=top_inc_unfair)) +
    geom_bar(aes(stat="identity")) +
    facet_wrap(~income, ncol = 5) + 
    scale_x_continuous(breaks=-4:4) +
    xlab("Perceived Unfairness of Top 10% Income") + ylab("Count"))

ggsave("figure/FigureB2_Appendix.jpeg", b2, width = 10, height = 5)

# Figure B3
(b3 <- dat %>%
    filter(is.na(income) == 0, is.na(job_self_unfair) == 0) %>%
    ggplot(aes(x=job_self_unfair)) +
    geom_bar(aes(stat="identity")) +
    facet_wrap(~income, ncol = 5) +
    scale_x_continuous(breaks=0:10) +
    xlab("Perceived Unfairness of One's Own Job Opportunities") + ylab("Count"))

ggsave("figure/FigureB3_Appendix.jpeg", b3, width = 10, height = 5)

# Figure A4
(b4 <- dat %>%
    filter(is.na(income) == 0, is.na(job_other_unfair) == 0) %>%
    ggplot(aes(x=job_other_unfair)) +
    geom_bar(aes(stat="identity")) +
    facet_wrap(~income, ncol = 5) +
    scale_x_continuous(breaks=0:10) +
    xlab("Perceived Unfairness of Others' Job Opportunities") + ylab("Count"))

ggsave("figure/FigureB4_Appendix.jpeg", b4, width = 10, height = 5)

##########################
# Figure 2
##########################
f2_manu_pred <- ggpredict(t1_1, c("netjobreduction_manu_scaled"))
(f2_manu <- f2_manu_pred %>%
    ggplot(aes(x=x, y=predicted)) +
    geom_line() +
    geom_ribbon(aes(ymin=conf.low, ymax=conf.high), alpha = 0.2) +
    theme(legend.position = c(0.8, 0.8)) + 
    scale_x_continuous(breaks = -8:8) +
    xlab("Net Reduction of Local Manufacturing Jobs") +
    ylab("Unfairness of Income for Self") )

f2_manu_job_pred <- ggpredict(t1_4, c("netjobreduction_manu_scaled"))
(f2_manu_job <- f2_manu_job_pred %>%
    ggplot(aes(x=x, y=predicted)) +
    geom_line() +
    geom_ribbon(aes(ymin=conf.low, ymax=conf.high), alpha = 0.2) +
    theme(legend.position = c(0.8, 0.8)) + 
    scale_x_continuous(breaks = -8:8) +
    xlab("Net Reduction of Local Manufacturing Jobs") +
    ylab("Unfairness of Job Opportunities for Self") )

f2_rti_pred <- ggpredict(t1_4, c("RTI")) 
(f2_rti <- f2_rti_pred %>%
    ggplot(aes(x=x, y=predicted)) +
    geom_line() +
    geom_ribbon(aes(ymin=conf.low, ymax=conf.high), alpha = 0.2) +
    scale_x_continuous(breaks = -2:2.5) +
    xlab("Routine Task Intensity") +
    ylab("Unfairness of Job Opportunities for Self") )

(f2 <- ggarrange(f2_manu, f2_manu_job, f2_rti, ncol = 3))
ggsave("figure/Figure2.jpeg", f2, width = 15, height = 5)

# Function for plotting
pred_plot <- function(pred_dat, xlab, ylab, scale, label = NULL, scale.y = NULL, label.y = NULL, ylim = NULL,
                      fontsize = 18) {
  if (is.null(ylim)==1){
    p <- pred_dat %>% 
      ggplot(aes(x=x, y=predicted)) +
      geom_line() +
      geom_ribbon(aes(ymin=conf.low, ymax=conf.high), alpha = 0.2) +
      theme_bw(base_size = fontsize) + scale_x_continuous(breaks = scale, labels = label.x) + 
      theme(axis.line = element_line(color='black'),
            plot.background = element_blank(),
            panel.grid.minor = element_blank(),
            panel.grid.major = element_blank()) +
      xlab(xlab) + ylab(ylab)
    return(p) 
  } else {
    p <- pred_dat %>% 
      ggplot(aes(x=x, y=predicted)) +
      geom_line() +
      geom_ribbon(aes(ymin=conf.low, ymax=conf.high), alpha = 0.2) +
      theme_bw(base_size = fontsize) + scale_x_continuous(breaks = scale, labels = label) + 
      theme(axis.line = element_line(color='black'),
            plot.background = element_blank(),
            panel.grid.minor = element_blank(),
            panel.grid.major = element_blank()) +
      scale_y_continuous(breaks = scale.y, labels = label.y, limits = ylim)+
      #ylim(ylim) +
      xlab(xlab) + ylab(ylab)
    return(p) 
  }
}

##########################
# Figure 3
##########################
f3_1_pred <- ggpredict(t2_1, c("inc_unfair [0:4]")) 
f3_2_pred <- ggpredict(t2_1, c("top_inc_unfair [0:4]")) 
f3_3_pred <- ggpredict(t2_1, c("bttm_inc_unfair [0:4]")) 

(f3_1 <- pred_plot(pred_dat = f3_1_pred, 
                   xlab = "Unfairness of Income for Self", 
                   ylab = "Support for Redistribution",
                   scale = 0:4, label = 0:4,
                   scale.y = seq(3.4, 4.0, 0.1), label.y = seq(3.4, 4, 0.1), ylim = c(3.4, 4)))

(f3_2 <- pred_plot(pred_dat = f3_2_pred, 
                   xlab = "Unfairness of Top 10% Income", 
                   ylab = "Support for Redistribution", scale = 0:4, label = 0:4,
                   scale.y = seq(3.4, 4.0, 0.1), label.y = seq(3.4, 4, 0.1), ylim = c(3.4, 4)))

(f3_3 <- pred_plot(pred_dat = f3_3_pred,
                   xlab = "Unfairness of Bottom 10% Income", 
                   ylab = "Support for Redistribution", scale = 0:4, label = 0:4,
                   scale.y = seq(3.4, 4.0, 0.1), label.y = seq(3.4, 4, 0.1), ylim = c(3.4, 4)))

f3 <- ggarrange(f3_1,f3_2,f3_3, ncol=3)
ggsave("figure/Figure3.jpeg", f3, width = 15, height = 5)

##########################
# Figure 4
##########################
f4_1_pred <- ggpredict(t2_3, c("inc_unfair [0:4]")) 
f4_2_pred <- ggpredict(t2_3, c("top_inc_unfair [0:4]")) 
f4_3_pred <- ggpredict(t2_3, c("bttm_inc_unfair [0:4]")) 

(f4_1 <- pred_plot(pred_dat = f4_1_pred, 
                   xlab = "Unfairness of Income for Self", 
                   ylab = "Support for Immigration", scale = 0:4, label = 0:4, 
                   scale.y = seq(3.7, 4.9, 0.2), label.y = seq(3.7, 4.9, 0.2), 
                   ylim = c(3.7, 4.9)))

(f4_2 <- pred_plot(pred_dat = f4_2_pred, 
                   xlab = "Unfairness of Top 10% Income", 
                   ylab = "Support for Immigration", scale = 0:4, label = 0:4, 
                   scale.y = seq(3.7, 4.9, 0.2), label.y = seq(3.7, 4.9, 0.2), ylim = c(3.7, 4.9)))

(f4_3 <- pred_plot(pred_dat = f4_3_pred, 
                   xlab = "Unfairness of Bottom 10% Income", 
                   ylab = "Support for Immigration", scale = 0:4, label = 0:4,
                   scale.y = seq(3.7, 4.9, 0.2), label.y = seq(3.7, 4.9, 0.2), ylim = c(3.7, 4.9)))

f4 <- ggarrange(f4_1,f4_2,f4_3, ncol=3)
ggsave("figure/Figure4.jpeg", f4, width = 15, height = 5)

##########################
# Figure 5
##########################
f5_1_pred <- ggpredict(t3_1, c("inc_unfair [0:4]")) 
f5_2_pred <- ggpredict(t3_1, c("top_inc_unfair [0:4]")) 
f5_3_pred <- ggpredict(t3_1, c("bttm_inc_unfair [0:4]")) 

(f5_1 <- pred_plot(pred_dat = f5_1_pred, 
                   xlab = "Unfairness of Income for Self", 
                   ylab = "Prob of Vote for the Radical Left",
                   scale = 0:4, label = 0:4,
                   scale.y = seq(0.02, 0.08, 0.01), label.y = seq(0.02, 0.08, 0.01), ylim = c(0.025, 0.085)))

(f5_2 <- pred_plot(pred_dat = f5_2_pred, 
                   xlab = "Unfairness of Top 10% Income", 
                   ylab = "Prob of Vote for the Radical Left",
                   scale = 0:4, label = 0:4,
                   scale.y = seq(0.02, 0.08, 0.01), label.y = seq(0.02, 0.08, 0.01), ylim = c(0.025, 0.085)))

(f5_3 <- pred_plot(pred_dat = f5_3_pred, 
                   xlab = "Unfairness of Bottom 10% Income", 
                   ylab = "Prob of Vote for the Radical Left",
                   scale = 0:4, label = 0:4,
                   scale.y = seq(0.02, 0.08, 0.01), label.y = seq(0.02, 0.08, 0.01), ylim = c(0.025, 0.085)))

f5 <- ggarrange(f5_1,f5_2,f5_3, ncol=3)
ggsave("figure/Figure5.jpeg", f5, width = 15, height = 5)

##########################
# Figure 6
##########################
f6_1_pred <- ggpredict(t3_3, c("inc_unfair [0:4]")) 
f6_2_pred <- ggpredict(t3_3, c("top_inc_unfair [0:4]")) 
f6_3_pred <- ggpredict(t3_3, c("bttm_inc_unfair [0:4]")) 

(f6_1 <- pred_plot(pred_dat = f6_1_pred, 
                   xlab = "Unfairness of Income for Self", 
                   ylab = "Prob of Vote for the Radical Right",
                   scale = 0:4, label = 0:4,
                   scale.y = seq(0.05, 0.17, 0.02), label.y = seq(0.05, 0.17, 0.02), ylim = c(0.06, 0.17)))

(f6_2 <- pred_plot(pred_dat = f6_2_pred, 
                   xlab = "Unfairness of Top 10% Income", 
                   ylab = "Prob of Vote for the Radical Right",
                   scale = 0:4, label = 0:4,
                   scale.y = seq(0.05, 0.17, 0.02), label.y = seq(0.05, 0.17, 0.02), ylim = c(0.06, 0.17)))

(f6_3 <- pred_plot(pred_dat = f6_3_pred, 
                   xlab = "Unfairness of Bottom 10% Income", 
                   ylab = "Prob of Vote for the Radical Right",
                   scale = 0:4, label = 0:4,
                   scale.y = seq(0.05, 0.17, 0.02), label.y = seq(0.05, 0.17, 0.02), ylim = c(0.06, 0.17)))

f6 <- ggarrange(f6_1,f6_2,f6_3, ncol=3)
ggsave("figure/Figure6.jpeg", f6, width = 15, height = 5)

##########################
# Figure 7
##########################
# Trust in Politics
f1_1 <- lm( trust_pols ~ inc_unfair + top_inc_unfair + bttm_inc_unfair + 
              college + female + age + income + union + unemployed + 
              citizen + religion + urban + cntry,
            data = dat,
            weights = dweight) 
f1_1 <- robustify(f1_1 , cluster = cntry)

f1_2 <- lm( sat_democracy ~ inc_unfair + top_inc_unfair + bttm_inc_unfair + 
              college + female + age + income + union + unemployed + 
              citizen + religion + urban + cntry,
            data = dat,
            weights = dweight) 
f1_2 <- robustify(f1_2 , cluster = cntry)

f1_3 <- lm( trust_pols ~ job_self_unfair + job_other_unfair +
              college + female + age + income + union + unemployed + 
              citizen + religion + urban + cntry,
            data = dat,
            weights = dweight) 
f1_3 <- robustify(f1_3, cluster = cntry)

# Panel 1 
f7_p1_self <- ggpredict(f1_1, c("inc_unfair [0:4]")) 
f7_p1_bottom <- ggpredict(f1_1, c("bttm_inc_unfair [0:4]")) 

f7_p1_bottom$group <- as.factor(2)
f7_p1 <- rbind(f7_p1_self, f7_p1_bottom)

(f7_1 <- f7_p1 %>%
  mutate(group = ifelse(group==1, "One's Own Income", "Bottom 10 %")) %>%
  ggplot(aes(x=x, y=predicted, linetype = group)) +
  geom_line() +
  geom_ribbon(aes(ymin=conf.low, ymax=conf.high), alpha = 0.2) +
  scale_y_continuous(limits = c(3,6.9)) + 
  theme_bw(base_size = 20) +
  theme(axis.line = element_line(color='black'),
        plot.background = element_blank(),
        panel.grid.minor = element_blank(),
        panel.grid.major = element_blank(),
        legend.position = c(0.35,0.15),
        legend.title = element_blank()) +
  xlab("Unfairness") + ylab("Trust in Politicians/Parliament"))

# Panel 2
f7_p2_self <- ggpredict(f1_2, c("inc_unfair [0:4]")) 
f7_p2_bottom <- ggpredict(f1_2, c("bttm_inc_unfair [0:4]")) 

f7_p2_bottom$group <- as.factor(2)
f7_p2 <- rbind(f7_p2_self, f7_p2_bottom)

(f7_2 <- f7_p2 %>%
  mutate(group = ifelse(group==1, "One's Own Income", "Bottom 10 %")) %>%
  ggplot(aes(x=x, y=predicted, linetype = group)) +
  geom_line() +
  geom_ribbon(aes(ymin=conf.low, ymax=conf.high), alpha = 0.2) +
  scale_y_continuous(limits = c(3,6.9)) + 
  theme_bw(base_size = 20) +
  theme(axis.line = element_line(color='black'),
        plot.background = element_blank(),
        panel.grid.minor = element_blank(),
        panel.grid.major = element_blank(),
        legend.position = c(0.35,0.15),
        legend.title = element_blank()) +
  xlab("Unfairness") + ylab("Satisfaction with Democracy"))

# Panel 3
f7_p3_self <- ggpredict(f1_3, c("job_self_unfair")) 
f7_p3_others <- ggpredict(f1_3, c("job_other_unfair")) 

f7_p3_others$group <- as.factor(2)
f7_p3 <- rbind(f7_p3_self, f7_p3_others)

(f7_3 <- f7_p3 %>%
  mutate(group = ifelse(group==1, "Job chance for self", "Job chance for others")) %>%
  ggplot(aes(x=x, y=predicted, linetype = group)) +
  geom_line() +
  geom_ribbon(aes(ymin=conf.low, ymax=conf.high), alpha = 0.2) +
  scale_y_continuous(limits = c(3,6.9)) + 
  theme_bw(base_size = 20) +
  theme(axis.line = element_line(color='black'),
        plot.background = element_blank(),
        panel.grid.minor = element_blank(),
        panel.grid.major = element_blank(),
        legend.position = c(0.4,0.83),
        legend.title = element_blank()) +
  xlab("Unfairness") + ylab("Trust in Politicians/Parliament"))

(f7 <- ggarrange(f7_1,f7_2,f7_3, ncol=3))

ggsave("figure/Figure7.jpeg", f7, width = 15, height = 5)

# New
(f7_1 <- f7_p1 %>%
    mutate(group = ifelse(group==1, "One's Own Income", "Bottom 10 %")) %>%
    ggplot(aes(x=x, y=predicted, linetype = group)) +
    geom_line() +
    geom_ribbon(aes(ymin=conf.low, ymax=conf.high), alpha = 0.2) +
    scale_y_continuous(limits = c(3.8,5.5)) + 
    theme_bw(base_size = 20) +
    theme(axis.line = element_line(color='black'),
          plot.background = element_blank(),
          panel.grid.minor = element_blank(),
          panel.grid.major = element_blank(),
          legend.position = c(0.35,0.15),
          legend.title = element_blank()) +
    xlab("Unfairness") + ylab("Trust in Politicians/Parliament"))

(f7_2 <- f7_p2 %>%
    mutate(group = ifelse(group==1, "One's Own Income", "Bottom 10 %")) %>%
    ggplot(aes(x=x, y=predicted, linetype = group)) +
    geom_line() +
    geom_ribbon(aes(ymin=conf.low, ymax=conf.high), alpha = 0.2) +
    scale_y_continuous(limits = c(5.4,7)) + 
    theme_bw(base_size = 20) +
    theme(axis.line = element_line(color='black'),
          plot.background = element_blank(),
          panel.grid.minor = element_blank(),
          panel.grid.major = element_blank(),
          legend.position = c(0.35,0.15),
          legend.title = element_blank()) +
    xlab("Unfairness") + ylab("Satisfaction with Democracy"))

(f7_3 <- f7_p3 %>%
    mutate(group = ifelse(group==1, "Job chance for self", "Job chance for others")) %>%
    ggplot(aes(x=x, y=predicted, linetype = group)) +
    geom_line() +
    geom_ribbon(aes(ymin=conf.low, ymax=conf.high), alpha = 0.2) +
    scale_y_continuous(limits = c(3,5.5)) + 
    theme_bw(base_size = 20) +
    theme(axis.line = element_line(color='black'),
          plot.background = element_blank(),
          panel.grid.minor = element_blank(),
          panel.grid.major = element_blank(),
          legend.position = c(0.37,0.15),
          legend.title = element_blank()) +
    xlab("Unfairness") + ylab("Trust in Politicians/Parliament"))

(f7 <- ggarrange(f7_1,f7_2,f7_3, ncol=3))

ggsave("figure/Figure7_new.jpeg", f7, width = 15, height = 5)

